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In the course of analyzing the dynamic behavior of mechanical systems 
subjected to random excitations, we investigate the associated Fokker-Planck 
equation. We also discuss the relationship between the characteristics of 
the random excitation and the nonlinear intensity coefficients governed by 
the ■physical properties of the system. This relationship leads to some 
simplified methods for solving the response probability density of certain 
nonlinear systems. We present general solutions to a class of multidimen- 
sional problems with desirable constraints. The random motion of a single- 
mode mechanical oscillator with various nonlinear stiffnesses and a charged 
particle moving in an electromagnetic field are examples. Cosine-power 
and sech-power distributions are found to be associated with the steady 
state response of a tangent stiffness system and a hyperbolic tangent stiff 
ness system, respectively. When the influencing magnetic vector potential 
M is irrotational, the stationary probability for the moving particle in the 
6-dimensional response phase-space is statistically independent. 

I. INTRODUCTION 

Although the theory of stochastic processes has found wide appli- 
cations in information and communication sciences for many years, 
only recent advances in rocket propulsion and aerospace industries 
have made random vibration problems subjects of growing importance 
in mechanical and civil engineering. These problems involve struc- 
tural responses due to random loadings and are in general nonlinear 
resulting from large motions. 1 Such nonlinear random transformation 
problems often encountered in practice are generally memory-de- 
pendent; that is, the equations of motion are described by nonlinear 
differential equations. 2-4 Under the Markov and Smoluchowski as- 
sumptions, it has been shown that the probability density function 
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of a large class of random processes satisfies equations of the Fokker- 
Planck (F-P) type. 5-0 Recently Pawula showed that generalized 
Fokker-Planck equations can be derived for many cases with both 
these assumptions removed. 7 Many interesting problems with their 
governing equations of the Fokker-Planck type have been investigated 
by various researchers. Rosenbluth, and others, studied the Fokker- 
Planck equation for the distribution function for gases with an inverse- 
square particle interaction force; 10 van Kampen used an Fokker- 
Planck equations to describe the thermal fluctuations in linear and 
nonlinear systems; 11 Ariaratnam found the steady state response distri- 
bution for a class of nonlinear two-mode mechanical oscillators by 
applying certain constraints to decouple the governing Fokker-Planck 
equation; 12 and Hempstead and Lax used Fourier transform tech- 
niques to eliminate the phase variable from the Fokker-Planck equa- 
tion in the polar coordinates for a rotating-wave Van der pol oscillator 
and found the phase and amplitude spectra. 13 

The Fokker-Planck equation, satisfied by the random response 
probability density function of a dynamic system, is a parabolic 
partial differential equation which generally is rather difficult to solve. 
Although approximate results may be obtained by using the per- 
turbation and equivalent linearization tehcniques (for which brief 
accounts are given in Appendixes A and B), the formal solution 
yielded by the appropriate Fokker-Planck equation is still the most 
sought one. 14-17 In this paper, we investigate the relationship between 
the random input characteristics and the intensity coefficients of the re- 
sponse process as governed by the physical properties of the system. 
Based on this relationship, we establish theorems concerning classes of 
potential-type and uncoupled-type solutions to the multidimensional 
stationary Fokker-Planck equations. Then, we present simple methods, 
based on these theorems, for solving such classes of random transfor- 
mations and we describe the required natural restraints which justify 
applying these methods on physical grounds. 

As examples we analyze two separate cases, a simple mechanical 
oscillator with various nonlinear spring resistances and a charged 
particle moving in an electromagnetic field, which we also consider to 
be subject to random excitation. 

II. MARKOV PROCESSES AND RANDOMLY EXCITED DYNAMIC SYSTEMS 

A stochastic process which has no aftereffect is called a Markov 
process. If y(t) = [yi(t), y 2 (y), ■•• , yjv(t)] is such a process (where 
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bold type indicates a vector) we can write 

p(7i I y 2 , • • • , vaO = p(y, I y 2 ), N > 2 (i) 

and 

p(ji ,y», ••• ,j*) = p(Yi I y 2 )p(y 2 1 y.) • • • p(y*-i I y*)p(y*) (2) 

where p( | ) and p( ) represent the conditional probability density and 
the joint probability density, respectively; and yi = y(ti), ■ • • ,y N = y(t N ). 

In examining the motion of a dynamic system under purely random 
disturbance, it is found that the phase point y(t 2 ) of the system at time 
t 2 depends only on the phase-point position y(£,) at the previous time 
t x . Therefore, the trajectory of the phase- point of a system under purely 
random disturbance is described by a Markov process y(t) = [yi(t), 
y»(t)> ' " > y*r(t)] in the phase space. The components y { , i=l, 3, 5, 
• • • , N — 1 represent the generalized coordinates of the system and 
components y i + 1 represent the first time derivatives of ?/,• . These N 
variables completely defined the dynamic state of a viboratory system 
in the iV-dimensional phase-space. 

If such a Markov process y(t) = [yi(t), y 2 (t), • • • , Pir(t)] satisfies the 
following conditions 

(i) the Smoluchowski equation 

p(ji I y 2 , At) = J dy 3 p(y, | y, , t + At)p(y a | y a , t) (3) 

holds for every y x , y 2 and y 3 defined in the 2V-dimensional phase space, 
(ii) the higher order intensity coefficients vanish, that is, 

K.(y) = lim -—■ <(y A , - y) 8 > = for s ^ 3 (4) 

A<->0 LA I 

and the first and second intensity coefficients 

K x (y) = lim j- (y,, A( - y.) = A t (j, t) (5) 

K 2 (y) = lim j- <(y,.. AI - y,.)(y,.. A( - y,)> = B it (j, t) (6) 

Al—O Ait 

exist, where the symbol ( ) indicates ensemble average,* 

(in) Ai , B ti are continuous and bounded, 

(iv) dAi/dyj exist for every y,- , and are continuous and bounded, 



* A Markov process which satisfies condition (ii) is said to be continuous. 



2034 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1969 

(v) dBij/dyt and d a B ij /dy i dy i exist for every y { and y f , and are 
continuous and bounded, and 

(vi) Bij(y, t) form a positive definite matrix, 

then the conditional probability p(y | y, t) of the process y(t) satisfies 
the Fokker-PIanck equation 

and the initial condition 

N 

P(Yo I y, t) = II 5(y< - y i0 ) as t -» 0, (8) 

where y is the initial state of y and S represents the Dirac delta function. 
Consider a class of n-degree-of-freedom nonlinear discrete dynamic 
systems whose motions are defined by the following system of differ- 
ential equations 

Xi + a,z,[l + e,D,(zi , x 2 , ■■ ■ , x n , x, , x 2 , • • • , £„)] 

+ b ( x,[l + HiSffa ,x 2 ,---,x nt x 1 ,---,x„)] = Vi(t) 

»-l,2, ...,n (9) 

where a,- and &,- are linear damping and stiffness coefficients, respec- 
tively; €,• and a*,- are nonlinear parameters; D { and S { are nonlinear 
functions; and rji(t) are excitations and are, in general, random processes. 
It is always convenient to transform the motion in n-dimensional 
generalized coordinates space into a 2n-dimensional phase-space, that 
is, let y { = Xi and y i+1 = x { . Equation (9) can be written in a system 
of 2n first order differential equations as 

Vi = fi(Vi , y 2 , ■ • • , y 2n ) + vS) i = 1, 2, ■ ■ • , 2n. (10) 

Assuming {^(t)) = and faCtOwtfa)) = Si/M^i _ *a) for constant 
Dy and applying equations (4) through (6) to equation (9) , we obtain 

Ai(y, t) = y i+l , 
A i+1 (j, t) = -a { y i+1 [l + e i D i (y l , ■ • • , y 2n )] 

- biViil + fuSiiVi , ••• , yjh,)]j 

£,.(y, = o, 
s,-,(y, = o, 
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and 



J + AI 



B<+t.t+i(y, = lim "Xt^ I 'h(r) o.t 



- a { y i+1 [l + e.D,] At - b,y t [l + M.S.] A* 

= lim -±- H J i?,(r,)i?.(T 2 ) d Tl dr 2 + 0(AO/ 

= S u . 

From the above results we notice that the coefficients A i+ i are de- 
termined by the specific nonlinear damping and nonlinear spring 
functions while coefficients fig depend only upon the statistical proper- 
ties of the random excitation functions. Therefore, we can conclude 
that for the random response components Vi(t) in equation (10) if 
the limits of their first and second increments 

lim i <Ay«0)> and lim j- (AyM AyM 

A<-»0 ^A t AI-0 /-it 

exist, then the probability density function of the many random func- 
tions y lf • • • , y 2n satisfies the Fokker-Planck equation (7) . 
Equation (7) can be written as 

,-_£•§£- (id 

where 

o*(j) = A t (j)p - h E £■ [B«(y)p] (12) 

describe the components of a probability current vector 
G = (Gi , G 2 , • • • , G 2n ). 

Because the general solution to the above multidimensional para- 
bolic partial differential equation with arbitrary boundary condition 
is difficult and impracticable to find, no such attempt is made in this 
study. However, with certain constraints on the properties of the 
system as well as on the characteristics of the random input, equation 
(11) may be reduced to simpler and more readily solvable forms. Such 
forms are considered in Section III. 
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III. STATIONARY SOLUTION OF JOINT PROBABILITY DISTRIBUTION 

In the limiting case of equation (11) when t — > oo the transition 
probability density p(y | y) tends to a stationary joint probability 
density p,,(y) — p„t(yi , V* , • • ■ , ViJ which is independent of the 
initial conditions and approaches a stationary (steady state) value as 
the time of passage is sufficiently large. Under this stationary situation, 
equation (11) becomes 

g*|M = (13) 

or 

V„,.-G = 0, »'- 1,2, ... ,2n. (14) 

Therefore G may be regarded as being incompressible and there are no 
sources or sinks in the region R. Notice that since rotational probability 
flows can occur even for cases of zero boundary conditions, that is, 
the probability current G(Gi , G 2 , • • ■ , G 2 „) satisfies 

(7,(30 =0, i = 1, 2, • • • , 2n (15) 

on the boundary of the region R under consideration, G may not vanish 
within R. In a special case, however, discussed in Section 3.1, the 
current vector G vanishes in the whole region R. 

Two situations which readily yield steady state solutions to equation 
(10) are investigated in Section 3.1. 

3. l Potential Distribution 

Under the assumption that the probability current vector G vanishes 
everywhere in R, a solution of the potential form exp (-U) can be 
constructed where U is a Liapunov type potential function of the system. 

If the equations of motion of a dynamic system satisfy the following 
conditions: 

(i) Gi = for all i, 

(ii) the matrix [B ti ] is not singular, 



»s:?*(?t , - u ') 



dy a i \i tyi 



and 
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(iv) there is no probability flow through the boundary of R, 
then the solution to the steady state Fokker-Planck equation (13) is 

p.,(y) = C exp {- f'"" """ y z\ZD ai ^-2ZD ai A i ~\dy a } 

(16) 
where [D al ] = [B i} }-\ that is, 

and C is the normalization constant determined by 

/ ••• / p»t dyt , ■•• , d2/ 2 „ = 1. (17) 

«'2n-fold •> 

Equation (16) can be easily verified by a direct substitution of p(y) = 
exp (— U) into equation (13) to solve for U. A special case of interest is 
when the matrix [£,-,] is isotropic, that is, when 

[£*,] = B(y)[M. (18) 

The solution for such a case, which can be obtained by substituting 
equation (18) into (16), is 



p.,(y) = Cexp(-C7) 

A 4 



(19) 



1 /• V < • " ' ■ V ■ n 2n 

C/ = l0g BW- 2 / a ,.... a .„ ? 



*<y) 

where C is the usual normalization factor. 



^,- 



3.2 Uncoupled Distribution 

There are cases when all generalized response variables of a system 
in the 2n phase-space coordinates are independent of one another. 
There are two approaches by which the solutions can be achieved. 

3.2.1 Forward Approach 

The governing Fokker-Planck equation may be reduced to a form 
in which the final probability distribution function clearly shows a 
statistical independent character. For this type of motion it is some- 
times possible to find appropriate partial operators which, when linearly 
operated on functions of the type gi(iji)p + hi(y*)(dp/dyi)t generate 
an equivalent nonlinear partial differential equation. If equation (13) 



V 
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can be put in the form 

£ U [gMp + Hv.) §£] - o (20) 

where the coefficients L, are arbitrary first order partial differential 
operators, and if there exists a p(y) satisfying each 

gMp + MvO ^ - 0, (t - 1, 2, • • • , 2»), 

then by Gray's uniqueness theorem such p(y) is the unique solution of 
equation (20). 18 Such a solution is 

,(y) = e ft ex P [- f jggrfx,] (2i) 

and C is the normalization factor. 

Notice that previous investigators such as Ariaratnam and Klein 
had their problems satisfying equation (20) and therefore obtained 
their solutions in a form similar to equation (21) . 

3.2.2 Backward Approach 

Sometimes it is more convenient to work backward. By this pro- 
cedure, the statistically independent property is assumed in order to 
derive the desirable solution for the Fokker-Planck equation under 
investigation. This method gives a close insight into the physical prop- 
erties of the system and enables the natural boundary conditions to be 
deduced. These deduced boundary conditions provide the necessary- 
constraints for randomly excited systems which have independent dis- 
tributions in their response variables. 

If equation (13) satisfies the conditions (22) 

(i) [B i} ] is a constant diagonal matrix for i, j = 1, 2, • • •, 2n, and, 
(u) the first order intensity coefficients Ai are functions of y t and 
Vi+i only, and there are no cross terms in A i} that is, 

A, = A^to) + A iti+1 (y i+1 ) 

A i+l = A i+Ui (yi) + A i+ui+l (y i+l )j 

for i = 1,3, ••• , 2n - 1; (23) 

then by setting p, t (y) = II<"i P*(tfO m equation (13), we obtain 

2 h B " pM dy] Vi{Vi) 



= E 



LpM {Ai ' f + Ai - i+ ° By< +P(y) B Vi . 
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Dividing the above equation by p(y), we obtain 

where the ' denotes d/dy t . 

The above equation can be reorganized as 

7Z[ L 2 p, p, J 

• -1,3, ••• -Pi Pi + 1 J 

A sufficient solution p(y) for the above equation requires that it satisfies 
the following relations 

B±±VZ _ Pi A _ 4/ = for » - i 2, • • • , 2n (24) 

2 p, p, 

and 

v' 1 »' 1 

(ii) — -j = — fii±i -j = E = constant 

Pi A i+lii p i+1 Ai, i+i 

for i - 1,3, ••• ,2n - 1. (25) 
Equation (24) can be reduced to 

Integration with respect to y { yields 

— pi = p,A,-,, +c, . 

When compared with equation (25) , c, vanishes and the following con- 
dition must hold 

A 2A, 'i = - A 2Ai n i+1 = E > (f - 1, 3, ■ • - , 2» - 1). (26) 

** i + 1 , iP i , i ** i . i + 1 D i + 1 , i + 1 

Equation (26) is the natural restraint under which the backward 
method can be applied. From the above analysis and the uniqueness 
theorem stated, we can claim. If the steady state Fokker-Planck 
equation satisfies conditions (22), (23), and (26), and if the phase-space 
coordinates y^ , y 2 , • • • , y 2n are statistically independent, then the 
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unique solution p.,(y) of equation (13) is 

f ' EA i+l , t (\i) dK 

-Jo 

- j"" +1 EA iti+1 (\ i+l ) d\ i+l j (27) 



2n-l 

P.i(y) - C II exp 

t' = 1.3, ••• 



where C is the normalized constant. 

rv. EXAMPLES 

4.1 Randomly Excited Nonlinear Simple Mechanical Oscillator 

When subjected to dynamic loadings rj(t) the equation of motion for 
a single-mode oscillator with nonlinear spring F(x) is 

x + fi± + F(a:) = v (t) (28) 

where the excitation q(t) is a gaussian stationary process with 

<*«> = j (29) 

in which S is a constant power spectral density. 

Letting y x = x and y 2 = x = fa , the intensity coefficients A { and 
Bi, can be found by using equations (4) and (5) as follows 

A, = A 1X + A 12 = y 2 , hence A u = 0, A 12 = 2/2 ] 

A 2 = A 2l + A 22 = -F{y x ) - Py 2 , 

hence A 21 = -/%,), A 22 = —0y 2 ; 

B 2 2 = 27T0 . 

Therefore the governing Fokker-Planck equation is 



dy 2 di/i di/ 2 

Forward approach: 

Equation 30 can be written in the form of equation (20) as 



(30) 



d \ „, , , rS dp 



_ dy 2 %i . 



?y 2 p + 



tS 6p 



= o, 



from which we see that 
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r d r d d 

ff 1(2/1) = ^(2/1), ^1(2/1) = ^ ' 



and 



£2(2/2) = 2/2 , ^2(2/2) = V s 



Substitution of g 1? <7 2 , h\ and /i 2 into equation (21) gives the following 
steady state solution for the system described by equation (28) 

*«-*-{-JT'9 l *-f£*} 

=^{-i[f™* + f} (3i) 

Backward approach: 

In this two-dimensional case, it can be shown that condition (26) is 
satisfied: 

E = Mii = ~2A 22 = -2(-gy 2 ) = _£_ # 
BuA 21 B 22 A l2 2irS a y 2 tS„ 

Therefore, according equations (22) through (27) the stationary solu- 
tion can be written : 

p..G/i , 2/ 2 ) = C exp f EA ai ({J #1 - f #4 12 (f 2 ) ^f 2 

-''O •'0 J 

[/; Fa,) *. + fl} 



= C exp i «■ 

7TO„ 



which is the same as equation (31) . 

4.1.1 Cubic Stiffness Characteristics 

Let us consider various nonlinear spring resistance functions F(x). 
The group of cubic stiffness characteristics is the classical case con- 
cerned with the hardening spring type of nonlinearity, which is repre- 
sented by 

F(x) = k x + ex 3 , (32) 

where k is the initial linear stiffness and e is the nonlinear coefficient. 
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Substituting equation (32) into (31), we obtain 



p, t (x, x) = C exp — 



L ttS \ 2 



.-„3T ex x 
2 + 4 + 2 



(33) 



Notice from equation (33) that the marginal probability density disitr- 
butions for x and x are statistically independent and that p(x) follows 
the gaussian law. As e — » 0, the system becomes linear and p(x) ap- 
proaches gaussian. 

4.1.2 Tangent Stiffness Characteristics™ 

The spring function F{x) in this case is shown in Fig. 1 and is repre- 
sented by 



_,, N (2k B d\ . 
F(x) = I — — ] tan 



wx 
2d) ' 



— d < x < d 



(34) 



in which m and d are constants. Notice that d is the threshold of the 
oscillator and m can be regarded as the mass of the oscillator. Again it 
follows from equation (31) that 



P. t (x, x) = p(x)p(x) 

- c [ exp (- *5z) 

F(I) 



4d 2 , irx 

exp \^2~2 m cos ^j 



(35) 



-d 




^J-'-'f'Kc 



Fig. 1 — Tangent stiffness characteristics. 
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where u] = kjm and a\ = wS„/20u>l is the corresponding linear mean 
square response [that is, if F(x) = k x]. 

The cosine-power distribution p{x) is shown in Fig. 2 for various 
values of n = 4d 2 /7r 2 o- 2 . It is interesting that for fixed u\ , p{x) ap- 
proaches the gaussian distribution as d — > °o and approaches the uniform 
distribution as d — > 0. 

4.1.3 Hyperbolic Tangent Stiffness Characteristics 

Figure 3 shows a hyperbolic tangent stiffness model; the spring 
resistance function F (x) is given by 

F(x) = -^rtanh (fix), b, k , m > 0. (36) 

mo 

Notice that the resistance F(x) developed during the vibration is 
bounded between kjbm and —k /bm. Therefore k /bm may be re- 
garded as the yield level and 1/bm the corresponding yielding response; 
equation (36) can be used to model the familiar elastic-perfect-plastic 
behavior observed in many physical realities. The joint probability 
density function for this can be obtained from equation (31), 



P.t(x, x) = p(x)p(x) = C exp - ^4^2 - pp In cosh (bx) \- 



(37) 



There are some interesting features about the limiting behavior of 
the marginal distribution p(x) which is of sech-power type 

p(x) = C 1 (6)sech (1/ff0,,,3> bx. (38) 

The sech-power term in equation (38) belongs to a class of distribution 
function f n (x) of a monotone decreasing sequence, integrable on [— °o , oo ] 
such that f n (x) ^ 1 for all x. It can be shown that 

lim j n (x) / J U(x) dx = 0; 20 and for lim p(x) = ( ^ exp [—^pj, 

the limit of the distribution function in equation (38) at zero, that is, 
lim t _ p(z) converges positively to a normal distribution with zero 
mean and variance <r 2 . The sech-power distribution p(x) is shown in 
Fig. 4 for various values of mf = l/<r z b 2 . 

4.2 Random Vibration of a Charged Particle Moving in an Electromag- 
netic Field 

As a second example, we consider a particle of mass m carrying 
charge q subjected to a random loading n(t) = rji(t)t + r) 2 {t)j + 77 3 (i)« 
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Fig. 2 — Cosine-power probability density distribution function. 

where i, j, & are base vectors in the Cartesian coordinates. The particle, 
moving in an electromagnetic field M = M y (x, y, z)i + M 2 (x, y, z)j + 
M a (x, y, z)fc, is also subjected to a friction force F, = — V.SF where 
^ = K A i y T ~f" a 2 w^ + X 3 y,) is an energy dissipation function. The com- 
plete force on the particle is therefore, 

F = gj- V<f>- \ ^ + \ (v X V X M)| (39) 

where <f> is the scalar potential and c is a constant. 

If we introduce a velocity-dependent potential 10, such that 

w = q<}> - 5 M-v (40) 

c 

then equation (39) can be written as 

F = -V d w + j t V,w. (41) 

Therefore the Lagrangian function L can be expressed in terms of w as 



L = T - q<f> + ^M-v 
c 



(42) 



where T represents the kinetic energy of the particle, and the equation 
of motion of the charged particle can be derived from the Lagrange 
equation 



A. ( dL \ _ Ok. _l i!L - n\ 
dt \dgj d qi + dg, " vA) 



(43) 
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Fig. 3 — Hyperbolic tangent stiffness characteristics. 



where g, represent the generalized coordinates of the motion. Using 
equations (42) and (43), one obtains 

a<£ . q (dM x . dM, 
dWi c \ dt du 3 

, dM t dM 2 dM 3 \ . . 

du 5 dUi dUi I 




-10 -a -6 -4 -2 



Fig. 4 — Scch-power probability density distribution function. 
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... 8<t> , q (dM 2 . 8M 2 

mu 3 + q— + * I — - + — - u 2 
8u 3 c \ dt 8u x 



(44) 



, 6M 2 dM 1 8M 3 



du* 



... d<f> , a (dM h . 



8u 3 
dM 3 



dth 



, dM 3 dM x dM 2 

+ -rr - u 2 — ^— w 2 — ^— w 4 ) + Mo = 7? a 



aw! 



aw. 



dUg 



where in equation (44), transformations u x = x, u 2 = x = Wj , u 3 = y, 
U4 = y = u 3 , u 5 = z, and u 6 = z = u 5 have been made. 
Define 



D t = (d 13 — d 2l )u 4 + (d 15 — d 31 )u 
D 2 = (d 21 — d 13 )u 2 + (d 25 — d 33 )u a 
D 3 = (d 31 — di 5 )u 2 + (d 33 — d 25 )u 4 J 



(45) 



where d aB = d aB (ui , u 3 , u 5 ) = 8M a /8u B , a = 1, 2, or 3 and /3 = 1, 3, 
or 5, are independent of the velocity variables. 

The following properties are assumed for the forcing function n(t) 



(v a ) = for a = 1, 2, or 3 

(vMVa(Q) = for a 9* p] 

= 2rS a (ti ~ it) for a = jSj 



(46) 
(47) 



From equations (44) through (47), the six-dimensional Fokker- 
Planck equation governing the transition probability density p(ux , 
u 2 , • • • , Ua), derived by the standard technique, is 



P(ui , ■ • • ,u a ) 

d ( \ d 



m m 8u Y mc \ dt 



d , . 8 j[ X 2 q 8<t> q J8M 2 . _ \"| \ 

~ 8V 3 M -^Al-m U '-mdV 3 -mc\-dT+ D2 )_ V l 



8u« 



dU a 



X, 



q 8<j> q(8M 3 



-•Z-M'-zr < -2«.-±?-i 



»I 



m dw ( 



mc 



\ d< 



+ 2), 






■ - a p , o d p „ d"p 
m L 8u 2 8u 4 8u 



] 



(48) 
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The initial condition for equation (48) is 

6 

p(«x ,u 2 , ■•• ,u e ) i,. = n •[■ 

i-l 

At t — > oo , the motion becomes stationary; th 
if dM/dt = 0. Imposing these conditions on 
use of equation (45), the corresponding s 
(48) is: 



refore p = and p = p„ 
Equation (48) and making 
;eady state equation for 



h. ±- _ I ±-)L UaV+ l[§L%L) +11 

m 2 du 2 m duj\ 2F \ y duj m dul 



, X 2 d 1 d \( tS 2 dp , 

\m du 4 m du 3 /\ r \ 2 duj 



(-, 



. X 3 <9 1 a V , ^ dp , 



.LJL f i£ tt^z dp 

n du t \ du 3 P X 2 du 3 

l_d_( d± ttS 3 dpi \ 
n dw„ V du 5 P X 3 duj 



+ 3 



(*£+*£+ ** M 



ap_ 



dU B 



= 0. 



Notice that in the above equation p = p, t ■ 
of equation (50) vanish if 



(50) 
Te :ms in the last parentheses 



(i) D, = D 2 = D 3 == 



or 



(w) w 2 : u 4 : w„ = (d 25 — d 33 ) : (d 31 — 



(in) V X M = or 



where G is a scalar potential function of x, y 

If any one of these conditions is satisfied 
same form as equation (20) and its solution 
by using equation (21), that is: 



p. t (u x , ■•■ , M 6 ) 



3f 



_ ah, f ' i£ d { 3 _ 

to 2 Jn dt 3 



mXa /w 



' • 5f 3 " 7rS 3 \2 

where C is the normalization factor determined by 



)LUTIONS 



(0 -!*<(«.)]. 
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(49) 



60 7T&! dp 

q dUi V X, dw, 



(51a) 



U : (d. 3 - d 2i ) (51b) 



M = G 



(51c) 



and 2. 

, equation (50) is of the 
lean be found immediately 



>-t(f) 



(52) 
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J ••■ J p„ du, , -■• , du a = 1. 



If clamping factors in three directions are identical and the random 
loading n is uniform, that is, if 

JF = \\{v\ + vl + i© (53) 

and 

0, = S 2 = S 3 = S, (54) 

then equation (52) becomes 

p., - Cexp[-^(T + #)]• (55) 

Applying this result to single-mode conservative oscillators with 
potential V (re) and subjected to gaussian white random loadings, the 
stationary response probability density is 

p., = Cexp(-^tf) (56) 

where H = T + V is the Hamiltonian function of the system. 

It is interesting that if the magnetic vector potential is irrotational, 
the steady state response probability density of a charged particle 
under white noise type random disturbances is statistically independ- 
ent. The solution can be immediately written down in terms of a 
quadrature. Extension of this result to a conservative dynamic sys- 
tem shows that stationary probability solution is of the form p Ht = 
C exp[— /(y, H)], where H is the Hamiltonian function of the system 
and y is a coefficient depending on the random input characteristics 
and the energy dissipation mechanism of the system. 

APPENDIX A 

Perturbation Technique 

The perturbation method is based on a series expansion in powers 
of the nonlinearity coefficient c. This method is valid for small values 
of c only. For example, consider a single-mode oscillator with non- 
linear damping and stiffness, when subjected to random force 7j(t) , the 
equation of motion is 

x + 2/34(1 + t £ b n x n ) + o&r.(l + * E anX n ) - v(t). (57) 



+ ?x n . 



FOKKER-PLANCK EQUATION S0LTJT; > 

A series solution can be assumed* 

x(t) = x„ + eXj + ex 2 + • • 
From equations (57) and (58) it follows that 

x + 2/3.r u + colx u = 77(0 

x, + 2/3i, + 0&:, = - 2 (2/3M" 

Therefore the correlation function is 

R xx (t) = B, ».(t) + efi IOIl (r) f d2,,..(T; 



where 

#*» - (x{t)y(t + r)j 

Equations (59) and (60) are linear and their 
obtained 



x 
and 



n (t) = f h(T) V (t - r) dr, 

Jo 



Xl (t) = -T h(r) Z [2fibX + \t - r) + 

•'0 n 



where h(r) = [e""7«.(l - ^ sin (1 - /3 2 )*co 
for system described by equation (59). 

From equations (63) and (64), the nonli 
be found. Considering only the first-order pefturbat 



R 



:....(*) = ~f Hr) £ [2flb n (x.(t)i n + \t - 

Jo n 

which can be evaluated explicitly if x {t) is 



APPENDIX B 

Equivalent Linearization Technique 
Consider the equation of motion 



x + 0x + u x + eg(x, x, t) 



♦Under fnirly general conditions it can be 
convergent. 
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(58) 



(59) 
1 +ulaX +1 ). (60) 



(61) 



(62) 
solutions can be readily 

(63) 



a n xT\t - r)] dr (64) 
,t is the transfer function 



in-, ar response moments can 
ion, 



r)> 

k(«K +1 (« " r))] dr, 

gaussian. 



= v(t), (65) 

showh that this series solution is 
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for a nonlinear function g, dependent on both x and x. The following 
equation is said to be equivalent to equation (65) in the sense that the 
mean square deficiency is minimized: 

x + 0.x + o a .x + e(x, x, t) = v (t). (66) 

The deficiency e(x, x, t) = (0 — (3 e )x + (col — w«)s + e 9( x > *i 0i m 
which 0, and u« are equivalent damping and frequency is determined by, 

-f£*- = and -£-§- = 0, (67) 

dp, dco, 

where the ( ) av indicates time average, that is, 

< e 2 ) av = limi f e\x,x,t)dt. (68) 

Using equations (67) and (68), we obtain 

2 2 | fo g/ay /7«\ 

w« — w t" e "7~2\ — v'UJ 

From equations (69) and (70) and by neglecting the deficiency term e, 
equation (66) can be solved by using standard linear theory. If the 
system is nonhereditary, the time averages in equations (67) through (70) 
are replaced by ensemble averages. /?. and co 2 , for this situation can be 
solved by a prior assumption for the probability density function of x(t). 
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